
setThreadOptions(numThreads=12,stackSize=1e07)
 
for(modeltype in c("simplified"#, "rich"
                   )){
  source("./auxfiles/loaddata.R")
  setwd(base_dir)
  
  
####PLOTS: EVALUATING FIT OF MODEL
######Targetted moments:
###TABLE OF ALL MOMENTS:

norun<-0
if(norun==0){
#setwd("./V_orig/")
   origtime<-Sys.time()
simdata<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                      shocks.wagereal,shocks.spwagereal,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                      marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=1)
newtime<-Sys.time()


saveRDS(simdata,"./modeloutput/simdata.RDS")

simdata<-cleansim(simdata,select=select,selectC=selectC)
#evaluating objective:
real_earnings<-data.table(read.dta13("permanentwage_byage.dta"))
evalfit(simdata=simdata,path="./modeloutput/",suffix=modeltype,onlyFT=1,nocomp=nocomp, nostockDI=nostockDI,noeventDI=noeventDI,nocontflow = nocontflow,
        shocks.wagenoise=shocks.wagenoise, shocks.spwagenoise=shocks.spwagenoise,
        consmoments_simple=moments_simple,
        workmoments_simple=moments_simple,
        DImoments_simple=moments_simple)
}

#full insurance condition:
#one QWC per year since quarter after turning 21.
#assuming you worked when 21 and 22 (ages not modeled).
simdata[,fullcredits:=cumsum(Work*4)+8,by=i]
simdata[,fullcredits_needed:=(age-21)/4]

#20/40 rule:
#assuming you worked when 21 and 22 (ages not modeled).
simdata[age<=31,agecredits:=cumsum(Work*4)+8,by=i]
for(a in 32:62){
  simdata[,agecreditsa:=sum(Work[age<=a & age>=a-10]),by=i]
  simdata[age==a,agecredits:=agecreditsa*4+8]
}
simdata[age<=23,agecredits_needed:=6]
simdata[age>=24 & age<=31,agecredits_needed:=(age - 21)*4/2]
simdata[age>31 ,agecredits_needed:=20]

simdata[Apply==1,mean(agecredits>=agecredits_needed & fullcredits>=fullcredits_needed)]
simdata[,meanW:=cumsum(Wage*(Work==1))/cumsum(Work),by=i]
simdata[yearDI==1,mean((meanW-Wage)/meanW)]

healthdata<-NULL
for(typeval in 1:3){
  ages <-0:49 + 23
  probs<-c(1,0,0,0,0,0)
  healthdata<-rbind(healthdata,
                    data.table(age = 23,
                               type=typeval,
                               H=t(probs)))
  for(t in 1:49){
    probs<-probs%*%healthprob[[typeval]][t,,]
    healthdata <- rbind(healthdata,
                        data.table(age=t+23,
                                   type=typeval,
                                   H=probs))
  }
}
healthdata[,sev:= H.V5 + H.V6]
healthdata[,mod:= H.V3 + H.V4]
healthdata[,any:= H.V3 + H.V4 + H.V5 + H.V6]
healthdata[,any:= H.V3 + H.V4 + H.V5 + H.V6]

pdf("./modeloutput/share_sevlimited.pdf")
print( ggplot(data=healthdata[age <= 62,])
       +geom_line(aes(x=age,y=sev,linetype=as.factor(type)),size=1)
       +scale_x_continuous(breaks=pretty_breaks(n=12))
       +scale_linetype_manual(values=c("solid","dashed","dotted"),labels=c("Low Type","Moderate Type","High Type"))
       +labs(x="Age",y="Share with Severe Work Limitation",linetype="")
       +theme(legend.position="bottom",
              legend.key.size=unit(1,'cm'))
       
)
dev.off()
pdf("./modeloutput/share_modlimited.pdf")
print( ggplot(data=healthdata[age <= 62,])
       +geom_line(aes(x=age,y=mod,linetype=as.factor(type)),size=1)
       +scale_x_continuous(breaks=pretty_breaks(n=12))
       +scale_linetype_manual(values=c("solid","dashed","dotted"),labels=c("Low Type","Moderate Type","High Type"))
       +labs(x="Age",y="Share with Moderate Work Limitation",linetype="")
       +theme(legend.position="bottom",
              legend.key.size=unit(1,'cm'))
       
)
dev.off()

pdf("./modeloutput/share_anylimited.pdf")
print( ggplot(data=healthdata[age <= 62,])
       +geom_line(aes(x=age,y=any,linetype=as.factor(type)),size=1)
       +scale_x_continuous(breaks=pretty_breaks(n=12))
       +scale_linetype_manual(values=c("solid","dashed","dotted"),labels=c("Low Type","Moderate Type","High Type"))
       +labs(x="Age",y="Share with Any Work Limitation",linetype="")
       +theme(legend.position="bottom",
              legend.key.size=unit(1,'cm'))
       
)
dev.off()
# rm(simdata_probs)
# gc()

marriagedata<-NULL
for(typeval in 1:3){
  ages <-0:49 + 23
  probs<-marriageprob_init[typeval]
  marriagedata<-rbind(marriagedata,
                      data.table(age = 23,
                                 type=typeval,
                                 M=probs))
  for(t in 1:49){
    probs<-probs*marriageprob[[typeval]][[1]][t,1] + (1-probs)*marriageprob[[typeval]][[1]][t,2]
    marriagedata <- rbind(marriagedata,
                          data.table(age=t+23,
                                     type=typeval,
                                     M=probs))
  }
}


pdf("./modeloutput/share_married.pdf")
print( ggplot(data=marriagedata[age <= 62,])
       +geom_line(aes(x=age,y=M,linetype=as.factor(type)),size=1)
       +scale_x_continuous(breaks=pretty_breaks(n=12))
       +scale_linetype_manual(values=c("solid","dashed","dotted"),labels=c("Low Type","Moderate Type","High Type"))
#       +ylim(c(0,1))
       +labs(x="Age",y="Share Married",linetype="")
       +theme(legend.position="bottom",
              legend.key.size=unit(1,'cm'))
       
)
dev.off()

shockdata<-read.dta13("./output/householdhealthshocks_stationary_marginal.dta")
shockdata<-as.data.table(shockdata)

pdf("./modeloutput/sevhealthshocks_stationary_marginal.pdf")
print( ggplot(data=shockdata[headsick==1,])
       +geom_line(aes(x=age,y=p,linetype=as.factor(type)),size=1)
       +scale_x_continuous(breaks=pretty_breaks(n=12))
       +scale_linetype_manual(values=c("solid","dashed","dotted"),labels=c("Low Type","Moderate Type","High Type"))
#       +ylim(c(0,1))
       +labs(x="Age",y="Prob. of Severe Health Shock for \n Healthy Ref. Person",linetype="")
       +theme(legend.position="bottom",
              legend.key.size=unit(1,'cm'))
)
dev.off()

###########EVENT STUDIES#############
#of disability onset, head:
#restricting to first observed onset of disability
#and matching to controls:
noevent<-0
if(noevent==0){
#  setwd("../eventstudy/figures")
#evalevent(simdata=simdata,path="figures/",onlyFT=1,eventtype="head")
evalevent(simdata=simdata,path="./modeloutput/",onlyFT=1,eventtype="headsev",suffix=modeltype)
#evalevent(simdata=simdata,path="figures/",onlyFT=1,eventtype="headmod")
}


#####

#####ELASTIFCITIES IMPLIED BY MODEL###
noelast<-0
if(noelast==0){
  shocks.wage<-shocks.wagereal
  shocks.spwage<-shocks.spwagereal
  
evalelast(path="./modeloutput/",paramvals=newguess, suffix=modeltype,
          marriageprob=marriageprob,
          marriageprob_init=marriageprob_init,
          onlyFT=1,workDI=0)
}

}

####

